pbmc6k clustering
Load and prepare the data
Load the the count matrix and remove all the genes/rows with zero counts.
# Load data.
pbmc.data <- read.table("~/Documents/Stage/Matrices/count-matrix-6k-v2.tsv", row.names=1, header=T)
# Remove all genes with zero counts.
pbmc.data <- pbmc.data[apply(pbmc.data[,-1], 1, function(x) !all(x==0)),]The rownames are ensembl gene ids, this is the standard output from featureCount. Convert the rownames from ensembl gene ids to gene name using biomaRt.
Get a dataframe with gene names corresponding to the ensemle id’s in the data.
library(biomaRt)
ensembl = useMart("ensembl", dataset="hsapiens_gene_ensembl")
gene.name.table <- getBM(attributes = c("external_gene_name", "ensembl_gene_id"),
filters = "ensembl_gene_id",
values = rownames(pbmc.data),
mart = ensembl)Not all ensembl ids are present in the BioMart database, so only set when available.
get.gene.name <- function(ensemble.id) {
gene.name <- gene.name.table[gene.name.table$ensembl_gene_id == ensemble.id, 1]
if (identical(gene.name, character(0)) || gene.name == "") {
gene.name <- ensemble.id
}
return(gene.name)
}
rownames(pbmc.data) <- make.names(lapply(rownames(pbmc.data), get.gene.name), unique = T)Searat can handle regular and sparse matrices, but the latter result in significant memory and speed savings. So create a sparse matrix.
pbmc.data.sparse <- Matrix(as.matrix(pbmc.data), sparse = T)Creating the searat object and QC
Initialize the Seurat object with the raw (non-normalized data). Keep all genes expressed in >= 3 cells, keep all cells with >= 200 genes. Perform log-normalization, first scaling each cell to a total of 1e4 molecules (as in Macosko et al. Cell 2015).
pbmc <- new("seurat", raw.data = pbmc.data.sparse)
pbmc <- Setup(pbmc, min.cells = 3, min.genes = 200, do.logNormalize = T, total.expr = 1e4, project = "10X_PBMC")Number of remaining cells = 5418/5419.
length(pbmc@cell.names)## [1] 5418Extract the mitochondrial genes.
mito.genes <- grep("^MT\\.", rownames(pbmc@data), value = T)
percent.mito <- colSums(expm1(pbmc@data[mito.genes, ]))/colSums(expm1(pbmc@data))
pbmc <- AddMetaData(pbmc, percent.mito, "percent.mito")VlnPlot(pbmc, c("nGene", "nUMI", "percent.mito"), nCol = 3)
Plot the total percentage of mitocondrial genes and the number of genes against the number of UMIs.
par(mfrow = c(1, 2))
GenePlot(pbmc, "nUMI", "percent.mito")
abline(h = 0.05, lty = 2)
GenePlot(pbmc, "nUMI", "nGene")
abline(h = 2500, lty = 2) 
par(mfrow = c(1, 1))Cut-off cells with cells with more then 2500 gene’s, they are possible multiplets. And cut-off cells with more than 0.05 percent mitochondrial rna.
pbmc <- SubsetData(pbmc, subset.name = "nGene", accept.high = 2500)
pbmc <- SubsetData(pbmc, subset.name = "percent.mito", accept.high = 0.05)
length(pbmc@cell.names)## [1] 51455142 cells remain.
pbmc@var.genes## logical(0)Regress out the number of molecules and the percentage mitochondrial genes.
pbmc <- RegressOut(pbmc, latent.vars = c("nUMI", "percent.mito"))Detect variable genes
Detection of variable genes across the single cells. This is done by “calculating the average expression and dispersion for each gene, placing these genes into bins, and then calculating a z-score for dispersion within each bin”. These are “typical parameter settings for UMI data that is normalized to a total of 1e4 molecules”.
pbmc <- MeanVarPlot(pbmc ,fxn.x = expMean, fxn.y = logVarDivMean, x.low.cutoff = 0.0125, x.high.cutoff = 3, y.cutoff = 0.5, do.contour = F)
length(pbmc@var.genes)PCA and tSNE
Only run PCA on variable genes. According to the toturial, using all the genes yields the same results, but is significantly slower.
pbmc <- PCA(pbmc, pc.genes = pbmc@var.genes, do.print = TRUE, pcs.print = 5, genes.print = 5)## [1] "PC1"
## [1] "TRBC2" "IL32" "CD3D" "IL7R" "CXCR4"
## [1] ""
## [1] "CST3" "FCN1" "TYROBP" "AIF1" "S100A8" "LST1"
## [1] ""
## [1] ""
## [1] "PC2"
## [1] "CD79A" "MS4A1" "HLA.DQA1" "CD79B" "IGHM"
## [1] ""
## [1] "NKG7" "CST7" "PRF1" "GZMA" "GZMB" "CTSW"
## [1] ""
## [1] ""
## [1] "PC3"
## [1] "GZMB" "GNLY" "FGFBP2" "PRF1" "NKG7"
## [1] ""
## [1] "IL7R" "CD3D" "VIM" "NOSIP" "MAL" "RGCC"
## [1] ""
## [1] ""
## [1] "PC4"
## [1] "CYBA" "TYROBP" "HLA.DPB1" "CTSS" "CD14"
## [1] ""
## [1] "PPBP" "PF4" "CLU" "TUBB1" "SDPR" "GNG11"
## [1] ""
## [1] ""
## [1] "PC5"
## [1] "FCGR3A" "MS4A7" "HES4" "CDKN1C" "CKB"
## [1] ""
## [1] "S100A8" "CD14" "GPX1" "MS4A6A" "S100A12" "VCAN"
## [1] ""
## [1] ""Project principal components on all genes (including non-variable genes).
pbmc <- ProjectPCA(pbmc)Examine and visualize PCA results a few different ways.
PrintPCA(pbmc, pcs.print = 1:5, genes.print = 10, use.full = TRUE)## [1] "PC1"
## [1] "RPS27A" "MALAT1" "PTPRCAP"
## [4] "RPSA" "RPL3" "RPS3A"
## [7] "ENSG00000229164" "RPS3" "TRBC2"
## [10] "IL32"
## [1] ""
## [1] "CST3" "LYZ" "S100A9" "FCN1" "TYROBP" "AIF1" "S100A8"
## [8] "LST1" "FTL" "FCER1G"
## [1] ""
## [1] ""
## [1] "PC2"
## [1] "CD79A" "MS4A1" "HLA.DQA1" "CD79B" "IGHM"
## [6] "IGHD" "TCL1A" "HLA.DRA" "LINC00926" "VPREB3"
## [1] ""
## [1] "NKG7" "CST7" "PRF1"
## [4] "GZMA" "GZMB" "CTSW"
## [7] "FGFBP2" "GNLY" "ENSG00000161570"
## [10] "GZMH"
## [1] ""
## [1] ""
## [1] "PC3"
## [1] "GZMB" "GNLY" "FGFBP2" "PRF1" "NKG7" "CST7"
## [7] "CD74" "HLA.DPB1" "GZMA" "SPON2"
## [1] ""
## [1] "ENSG00000229164" "IL7R" "RPS14"
## [4] "CD3D" "RPL32" "LDHB"
## [7] "RPS12" "RPL11" "VIM"
## [10] "JUNB"
## [1] ""
## [1] ""
## [1] "PC4"
## [1] "RPS2" "RPL13" "MALAT1" "RPL19" "RPS6" "RPL13A" "RPL10"
## [8] "RPL28" "RPS14" "RPL11"
## [1] ""
## [1] "PPBP" "PF4" "CLU" "TUBB1" "SDPR"
## [6] "GNG11" "ACRBP" "HIST1H2AC" "GP9" "TREML1"
## [1] ""
## [1] ""
## [1] "PC5"
## [1] "FCGR3A" "MS4A7" "HES4" "CDKN1C" "CKB"
## [6] "RHOC" "CTSL" "LINC01272" "MS4A4A" "SIGLEC10"
## [1] ""
## [1] "S100A8" "CD14" "GPX1" "MS4A6A" "S100A12" "VCAN" "LGALS2"
## [8] "FOLR3" "S100A9" "RNASE6"
## [1] ""
## [1] ""VizPCA(pbmc, 1:2)
PCAPlot(pbmc, 1, 2)
PCHeatmap(pbmc, pc.use = 1, cells.use = 100, do.balanced = TRUE)
PCHeatmap(pbmc, pc.use = 1:12, cells.use = 500, do.balanced = TRUE, label.columns = FALSE, use.full = FALSE)
To compute the significance of the PCs the tutorial uses the “JackStraw” method. This is very CPU intensice, I aborted after 1 hour. A more ad hoc method is to identify the elbow of this plot, and select the PC’s up to this points.
PCElbowPlot(pbmc)
The elbow falls around PC 10.
Determine distances bases on the first 10 PC’s. “resolutioan parameter between 0.6-1.2 typically returns good results for single cell dataset” save.SNN=T saves the SNN so that the SLM algorithm can be rerun using the same graph, but with a different resolution value.
pbmc <- FindClusters(pbmc, pc.use = 1:10, resolution = 0.6, print.output = 0, save.SNN = T)Run tSNE
pbmc <- RunTSNE(pbmc, dims.use = 1:10, do.fast = T)TSNEPlot(pbmc, do.label = T)
Identify marker genes
B-cells
b.cell.markers <- FindMarkers(pbmc, ident.1 = 3, min.pct = 0.25)
head(b.cell.markers, 20)## p_val avg_diff pct.1 pct.2
## LYZ 0 -3.263452 0.256 0.510
## CD79A 0 3.240544 0.927 0.024
## S100A9 0 -3.235467 0.108 0.366
## IGLC2 0 2.798172 0.368 0.017
## TCL1A 0 2.717464 0.564 0.012
## IGHM 0 2.631317 0.718 0.028
## CD79B 0 2.629724 0.911 0.117
## CST3 0 -2.603591 0.082 0.377
## NKG7 0 -2.601726 0.051 0.274
## MS4A1 0 2.574706 0.756 0.026
## TYROBP 0 -2.499070 0.076 0.412
## ENSG00000161570 0 -2.474276 0.083 0.284
## IGLC3 0 2.440580 0.306 0.012
## S100A4 0 -2.285672 0.285 0.855
## IGHD 0 2.193298 0.601 0.010
## S100A8 0 -2.178336 0.056 0.295
## IL32 0 -2.167829 0.058 0.579
## AIF1 0 -2.162039 0.057 0.399
## HLA.DQA1 0 2.161898 0.887 0.110
## CD74 0 2.157761 1.000 0.688Tc-cells
t.c.cell.markers <- FindMarkers(pbmc, ident.1 = c(5,6,8), min.pct = 0.25)
head(t.c.cell.markers, 20)## p_val avg_diff pct.1 pct.2
## LYZ 0 -3.638742 0.292 0.514
## S100A9 0 -3.299076 0.113 0.377
## CST3 0 -2.711585 0.119 0.383
## HLA.DRA 0 -2.662665 0.230 0.550
## TYROBP 0 -2.405417 0.094 0.424
## ENSG00000161570 0 2.354885 0.772 0.151
## S100A8 0 -2.222531 0.053 0.306
## FCER1G 0 -2.137002 0.062 0.402
## FCN1 0 -2.112921 0.059 0.323
## GZMK 0 1.989788 0.365 0.034
## AIF1 0 -1.914393 0.161 0.394
## GZMH 0 1.834503 0.325 0.056
## CD8A 0 1.789937 0.497 0.019
## CD8B 0 1.762056 0.533 0.028
## LST1 0 -1.747490 0.097 0.377
## CD74 0 -1.685492 0.655 0.745
## FTL 0 -1.587883 0.975 0.977
## TRGC2 0 1.560373 0.360 0.026
## CTSS 0 -1.474518 0.254 0.528
## FTH1 0 -1.460624 0.946 0.979Th-cells
t.h.cell.markers <- FindMarkers(pbmc, ident.1 = c(2,0,9), min.pct = 0.25)
head(t.h.cell.markers, 20)## p_val avg_diff pct.1 pct.2
## LYZ 0 -3.748129 0.304 0.577
## HLA.DRA 0 -3.460662 0.160 0.692
## S100A9 0 -3.457838 0.122 0.454
## CST3 0 -3.109594 0.081 0.488
## NKG7 0 -2.886006 0.059 0.353
## TYROBP 0 -2.861925 0.067 0.543
## CD74 0 -2.608642 0.516 0.855
## FCER1G 0 -2.501798 0.061 0.509
## S100A8 0 -2.474672 0.055 0.385
## FCN1 0 -2.402159 0.061 0.405
## HLA.DPA1 0 -2.357872 0.130 0.655
## HLA.DPB1 0 -2.309198 0.166 0.692
## ENSG00000161570 0 -2.186210 0.111 0.344
## HLA.DRB1 0 -2.158267 0.086 0.639
## AIF1 0 -2.137789 0.159 0.468
## LST1 0 -1.977795 0.107 0.460
## FTL 0 -1.930979 0.962 0.985
## LGALS2 0 -1.922605 0.023 0.318
## CTSS 0 -1.721246 0.259 0.611
## CFD 0 -1.598963 0.025 0.335Monocytes
monocyte.markers <- FindMarkers(pbmc, ident.1 = c(1,4), min.pct = 0.25)
head(monocyte.markers, 20)## p_val avg_diff pct.1 pct.2
## S100A9 0 3.964127 0.943 0.134
## LYZ 0 3.433368 0.994 0.309
## S100A8 0 3.230697 0.868 0.068
## FCN1 0 3.036080 0.908 0.075
## AIF1 0 2.866404 0.972 0.154
## CST3 0 2.815163 0.996 0.126
## TYROBP 0 2.704662 0.998 0.164
## LST1 0 2.620823 0.973 0.122
## ENSG00000161570 0 -2.579455 0.077 0.316
## FTL 0 2.531357 1.000 0.969
## ENSG00000229164 0 -2.516875 0.077 0.685
## LGALS2 0 2.462783 0.726 0.042
## NKG7 0 -2.349143 0.148 0.276
## CFD 0 2.329807 0.811 0.030
## LTB 0 -2.326188 0.262 0.786
## PTPRCAP 0 -2.287738 0.081 0.788
## FCER1G 0 2.286345 0.959 0.145
## IL32 0 -2.279541 0.071 0.651
## TRBC2 0 -2.225387 0.080 0.722
## IFITM3 0 2.187752 0.687 0.056Natural killer cells
nk.cell.markers <- FindMarkers(pbmc, ident.1 = 7, min.pct = 0.25)
head(nk.cell.markers, 20)## p_val avg_diff pct.1 pct.2
## GNLY 0 3.965651 0.977 0.082
## GZMB 0 3.289642 0.930 0.067
## S100A9 0 -2.922883 0.103 0.346
## LYZ 0 -2.894234 0.309 0.486
## FGFBP2 0 2.843483 0.834 0.056
## NKG7 0 2.735786 0.997 0.198
## PRF1 0 2.734806 0.927 0.091
## SPON2 0 2.353857 0.698 0.030
## HLA.DRA 0 -2.328911 0.203 0.513
## GZMA 0 2.243235 0.904 0.130
## CTSW 0 2.222552 0.957 0.212
## CST7 0 2.156877 0.930 0.139
## CST3 0 -2.094161 0.203 0.346
## HOPX 0 2.051723 0.645 0.091
## TRDC 0 2.025169 0.595 0.016
## ENSG00000129277 0 2.016787 0.615 0.065
## IGFBP7 0 1.969895 0.525 0.061
## S100A8 0 -1.961566 0.070 0.275
## FCN1 0 -1.923119 0.076 0.290
## CLIC3 0 1.877138 0.548 0.025Dendritic cells
dendritic.cell.markers <- FindMarkers(pbmc, ident.1 = 10, min.pct = 0.25)
head(dendritic.cell.markers, 20)## p_val avg_diff pct.1 pct.2
## FCER1A 0 2.869373 0.837 0.005
## ENSG00000229164 0 -2.113909 0.128 0.544
## CLEC10A 0 2.074088 0.686 0.011
## CD3D 0 -2.009190 0.047 0.490
## IL32 0 -1.913357 0.081 0.517
## HLA.DPB1 0 1.885365 1.000 0.489
## HLA.DPA1 0 1.833409 1.000 0.451
## TRBC2 0 -1.784020 0.186 0.572
## CST3 0 1.713351 0.988 0.327
## HLA.DRA 0 1.669098 1.000 0.486
## HLA.DQA1 0 1.646905 0.988 0.200
## HLA.DRB1 0 1.644483 1.000 0.425
## LTB 0 -1.564740 0.442 0.662
## CD3E 0 -1.562230 0.128 0.443
## CD74 0 1.505545 1.000 0.725
## HLA.DRB5 0 1.426911 0.988 0.320
## CD1C 0 1.392584 0.488 0.021
## NAPSB 0 1.312967 0.767 0.084
## PLD4 0 1.286553 0.488 0.033
## HLA.DRB6 0 1.270707 0.942 0.251Megakaryocytes
megakaryo.markers <- FindMarkers(pbmc, ident.1 = 11, min.pct = 0.25)
head(megakaryo.markers, 20)## p_val avg_diff pct.1 pct.2
## PPBP 0 5.773551 1.000 0.016
## PF4 0 4.392539 1.000 0.004
## SDPR 0 4.039698 0.871 0.006
## GNG11 0 4.035747 0.903 0.005
## CLU 0 3.626882 1.000 0.008
## HIST1H2AC 0 3.611592 0.935 0.022
## TUBB1 0 3.566211 0.871 0.006
## TSC22D1 0 3.440437 0.645 0.012
## ACRBP 0 3.141543 0.839 0.004
## PTCRA 0 3.125679 0.645 0.002
## RUFY1 0 3.084786 0.871 0.042
## TREML1 0 3.062469 0.742 0.001
## SPARC 0 3.059459 0.774 0.005
## ENSG00000138293 0 3.034532 0.839 0.121
## NRGN 0 2.993182 0.935 0.024
## CD9 0 2.982024 0.742 0.021
## GP9 0 2.981799 0.774 0.002
## CLEC1B 0 2.972835 0.581 0.001
## SNCA 0 2.779238 0.613 0.012
## MYL9 0 2.770657 0.645 0.002MS4A1 and CD79A are known markers for B-cells.
VlnPlot(pbmc, c("MS4A1", "CD79A")) 
Cluster 3 seems to be B-cells.
VlnPlot(pbmc, c("NKG7"), use.raw = T, y.log = T)
NKG7 is a marker for Natural killer cells but is expressed in T-cells also.
A few known marker genes
FeaturePlot(pbmc, c("MS4A1", "GNLY","CD3E","CD14","FCER1A","FCGR3A", "LYZ", "PPBP", "CD8A"), cols.use = c("grey","blue"))
Heatmaps can also be a good way to examine heterogeneity within/between clusters. The DoHeatmap() function will generate a heatmap for given cells and genes. In this case, we are plotting the top 20 markers (or all markers if less than 20) for each cluster.
pbmc.markers %>% group_by(cluster) %>% top_n(10, avg_diff) -> top10
DoHeatmap(pbmc, genes.use = top10$gene, order.by.ident = TRUE, slim.col.label = TRUE, remove.key = TRUE)
Cluster 0 and 1 look similar
Find the discriminating marker genes, between these clusters
cell.0.1.markers <- FindMarkers(pbmc, 0, 1)
head(cell.0.1.markers, 10)## p_val avg_diff pct.1 pct.2
## LYZ 0 -4.994542 0.305 1.000
## S100A9 0 -4.752301 0.137 0.986
## CST3 0 -4.002295 0.087 0.995
## TYROBP 0 -3.765259 0.073 0.997
## S100A8 0 -3.693892 0.054 0.987
## FCN1 0 -3.431628 0.060 0.969
## HLA.DRA 0 -3.165587 0.184 0.937
## LTB 0 3.067149 0.947 0.152
## AIF1 0 -3.066694 0.092 0.962
## FCER1G 0 -3.041993 0.057 0.945FeaturePlot(pbmc, c("ANXA2", "S100A4", "CCR7"), cols.use = c("green", "blue"))
cd3.genes <- c("CD3E", "CD3D", "CD3G")NK cell markers wikipedia: CD16+, CD56+ (NCAM1), CD3-, CD31, CD30, CD38 abcam: CD3-, CD56+ (NCAM1), CD94+, NKp46+ seurat: GNLY (Nk & T-cells), NKG7 (Nk & T-cells)
FeaturePlot(pbmc, c(cd3.genes, "NCAM1", "NKG7", "GNLY"), cols.use = c("green", "blue"), reduction.use = "tsne")
CD3 genes should be absent. NCAM1, NKG7 & GNLY are known markers
B cell markers CD45+(PTPRC) (also in T cells), CD19+, CD20+ (MS4A1), CD24+, CD38, CD22
grep("^CD22", rownames(pbmc@data), value = T)## [1] "CD226" "CD22"FeaturePlot(pbmc, c("CD19", "PTPRC", "MS4A1", "CD38", "CD22"), cols.use = c("green", "blue"), reduction.use = "tsne")
T cells CD3 and PTPRC/CD45(also B cells) are known markers van T cells
FeaturePlot(pbmc, c(cd3.genes, "PTPRC"), cols.use = c("green", "blue"), reduction.use = "tsne")
CD4+ T(helper) markers: CD4+ IL7R+ CD8+ T(cytotoxic) markers: CD8A, CD8B
FeaturePlot(pbmc, c("CD8A", "CD8B"), cols.use = c("green", "blue"), reduction.use = "tsne")
Monocytes
FeaturePlot(pbmc, c("CD14", "FCGR3B", "FCGR3A"), cols.use = c("green", "blue"), reduction.use = "tsne")
Macrophage CD11b+ (ITGAM), CD68+, CD163+
FeaturePlot(pbmc, c("ITGAM", "CD68", "CD163"), cols.use = c("green", "blue"), reduction.use = "tsne")
Build a network with B cells
# load("./seurat_pbmc_6k.RData")
#
# length(pbmc@data[,1])
#
# b.cells <- pbmc@scale.data[,WhichCells(pbmc, 3)]
# t.c.cells <- pbmc@scale.data[,WhichCells(pbmc, c(4,7))]
# t.h.cells <- pbmc@scale.data[,WhichCells(pbmc, c(0,1))]
# monocytes <- pbmc@scale.data[,WhichCells(pbmc, c(2,5))]
# nk.cells <- pbmc@scale.data[,WhichCells(pbmc, 6)]
# dendritic.cells <- pbmc@scale.data[,WhichCells(pbmc, 8)]
# megakaryos <- pbmc@scale.data[,WhichCells(pbmc, 9)]
#
# # CCR5 CTLA4 FOXP3 HLA-DQA1 HLA.DQB1 HLA.DRB1 HNF1A IL2RA IL6 INS ITPR3 OAS1 PTPN22 SUMO4
# # CCR5: is expressed by T cells and macrophages, and is known to be an important co-receptor for macrophage-tropic virus, including HIV, to enter host cells.
# # Certain HLA haplotypes are associated with a higher risk of developing type 1 diabetes, with particular combinations of HLA-DQA1, HLA-DQB1, and HLA-DRB1 gene variations resulting in the highest risk
# # CTLA4: Mutations in this gene have been associated with insulin-dependent diabetes mellitus.
#
# b.cells.cor <- cor(t(b.cells), method="pearson")
# t.c.cells.cor <- cor(t(t.c.cells), method="pearson")
# t.h.cells.cor <- cor(t(t.h.cells), method="pearson")
# monocytes.cor <- cor(t(monocytes), method="pearson")
# nk.cells.cor <- cor(t(nk.cells), method="pearson")
# dendritic.cells.cor <- cor(t(dendritic.cells), method="pearson")
# megakaryos.cor <- cor(t(megakaryos), method="pearson")
#
# diabetes.genes <- c("CCR5", "CTLA4", "FOXP3", "HLA.DQA1", "HLA.DQB1", "HLA.DRB1", "IL2RA", "IL6", "ITPR3", "OAS1", "PTPN22", "SUMO4")
#
# b.cells.diabetes <- pbmc@scale.data[diabetes.genes, WhichCells(pbmc, 3)]
# b.cells.dia.cor <- cor(t(b.cells.diabetes), method="pearson")
# b.cells.dia.cor * b.cells.dia.cor
# b.cells.dia.network <- GeneNetwork(data = b.cells.dia.cor, start.gene = "CCR5", min.correlation = 0.85)
# b.dia.net = network(b.cells.dia.network@network, directed = FALSE)
# ggnet2(b.dia.net, label = T) + ggtitle("B cells")
#
#
# source(file = "GeneNetwork.R")
#
# gene.symbol <- "CCR5"
# min.correlation <- 0.85
#
# #is.correlated <- (t.c.cells.cor[gene.symbol,] > min.correlation | t.c.cells.cor[gene.symbol,] < -min.correlation)
# is.correlated <- order(abs(b.cells.cor[gene.symbol,]),decreasing = T)[1:5]
# corelated.genes <- b.cells.cor[gene.symbol, is.correlated]
# corelated.genes
#
# b.cell.network <- GeneNetwork(data = b.cells.cor, start.gene = "CCR5", min.correlation = 0.85)
# b.cell.network <- addGene(b.cell.network, gene.name = "HLA.DQB1", min.correlation = 0.85)
# t.c.cell.network <- GeneNetwork(data = t.c.cells.cor, start.gene = "CCR5", min.correlation = 0.85)
# t.c.cell.network <- addGene(t.c.cell.network, gene.name = "HLA.DQB1", min.correlation = 0.85)
# t.h.cell.network <- GeneNetwork(data = t.h.cells.cor, start.gene = "CCR5", min.correlation = 0.85)
# t.h.cell.network <- addGene(t.h.cell.network, gene.name = "HLA.DQB1", min.correlation = 0.85)
# monocytes.network <- GeneNetwork(data = monocytes.cor, start.gene = "CCR5", min.correlation = 0.85)
# monocytes.network <- addGene(monocytes.network, gene.name = "HLA.DQB1", min.correlation = 0.85)
# nk.cells.network <- GeneNetwork(data = nk.cells.cor, start.gene = "CCR5", min.correlation = 0.85)
# nk.cells.network <- addGene(nk.cells.network, gene.name = "HLA.DQB1", min.correlation = 0.85)
# dendritic.cells.network <- GeneNetwork(data = dendritic.cells.cor, start.gene = "CCR5", min.correlation = 0.85)
# dendritic.cells.network <- addGene(dendritic.cells.network, gene.name = "HLA.DQB1", min.correlation = 0.85)
# megakaryos.network <- GeneNetwork(data = megakaryos.cor, start.gene = "CCR5", min.correlation = 0.85)
# megakaryos.network <- addGene(megakaryos.network, gene.name = "HLA.DQB1", min.correlation = 0.85)
#
# library(GGally)
# library(network)
# library(sna)
# library(ggplot2)
# require(gridExtra)
#
# b.net = network(b.cell.network@network, directed = FALSE)
# t.c.net = network(t.c.cell.network@network, directed = FALSE)
# t.h.net = network(t.h.cell.network@network, directed = FALSE)
# monocytes.net = network(monocytes.network@network, directed = FALSE)
# n.k.net = network(nk.cells.network@network, directed = FALSE)
# dendritic.net = network(dendritic.cells.network@network, directed = FALSE)
# megakaryos.net = network(megakaryos.network@network, directed = FALSE)
#
# b.plot <- ggnet2(b.net, label = T) + ggtitle("B cells")
# t.c.plot <- ggnet2(t.c.net, label = T) + ggtitle("T-c cells")
# t.h.plot <- ggnet2(t.h.net, label = T) + ggtitle("T-h cells")
# monocytes.plot <- ggnet2(monocytes.net, label = T) + ggtitle("Monocytes")
# n.k.plot <- ggnet2(n.k.net, label = T) + ggtitle("Natural killer cells")
# dendritic.plot <- ggnet2(dendritic.net, label = T) + ggtitle("Dentritic cells")
# megakaryos.plot <- ggnet2(megakaryos.net, label = T) + ggtitle("Megakaryocytes")
#
# b.plot;t.c.plot;t.h.plot
# monocytes.plot;n.k.plot;dendritic.plot;megakaryos.plot
#
# grid.arrange(b.plot, t.c.plot, t.h.plot, nrow=2, ncol=2)